Food competition analysis

Load packages

home <- here::here()

library(tidyverse)
library(tidylog)
library(RCurl)
library(sdmTMB)
library(RColorBrewer)
library(devtools)
library(patchwork)
library(ggstats)
library(ggh4x)
library(viridis)
library(sdmTMBextra)

# Source map-plot
#source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
source(paste0(home, "/R/functions/map-plot.R"))

Reading layer `StatRec_map_Areas_Full_20170124' from data source 
  `/Users/maxlindmark/Dropbox/Max work/R/cod-interactions/data/shapefiles/ICES-StatRec-mapto-ICES-Areas/StatRec_map_Areas_Full_20170124.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 11074 features and 17 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -44 ymin: 36 xmax: 69 ymax: 85
Geodetic CRS:  WGS 84

Read data & prepare data

d <- read_csv(paste0(home, "/data/clean/aggregated_stomach_data.csv")) |> 
  drop_na(group) |> 
  drop_na(oxy)

t <- d |> filter(fle_kg_km2 > 0)
t <- d |> filter(mcod_kg_km2 > 0)
t <- d |> filter(scod_kg_km2 > 0)

# Calculate relative prey weights (saduria and benthos)
d <- d |> 
  rename(oxygen = oxy) %>% 
  mutate(tot_weight = rowSums(select(., ends_with('_tot'))),  
         benthic_weight = amphipoda_tot + bivalvia_tot + gadus_morhua_tot +
           gobiidae_tot + mysidae_tot + non_bio_tot + 
           other_crustacea_tot + other_tot + other_pisces_tot + platichthys_flesus_tot +
           polychaeta_tot + saduria_entomon_tot) |> 
  rename(saduria_weight = saduria_entomon_tot,
         flounder_density = fle_kg_km2,
         large_cod_density = mcod_kg_km2,
         small_cod_density = scod_kg_km2) |> 
  mutate(tot_rel_weight = tot_weight / (pred_weight_g - tot_weight), 
         benthic_rel_weight = benthic_weight / (pred_weight_g - tot_weight),
         saduria_rel_weight = saduria_weight / (pred_weight_g - tot_weight)) |> 
  dplyr::select(-ends_with("_tot")) |> 
  dplyr::select(-predator_latin_name, date) |> 
  # Add small constant to large cod density because we want to take the log of it
  mutate(large_cod_density = ifelse(large_cod_density == 0,
                                    min(filter(d, mcod_kg_km2 > 0)$mcod_kg_km2)*0.5,
                                    large_cod_density),
         flounder_density = ifelse(flounder_density == 0,
                                   min(filter(d, fle_kg_km2 > 0)$fle_kg_km2)*0.5,
                                   flounder_density),
         small_cod_density = ifelse(small_cod_density == 0,
                                    min(filter(d, scod_kg_km2 > 0)$scod_kg_km2)*0.5,
                                    small_cod_density)) |> 
  # Scale variables
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         fhaul_id = as.factor(haul_id),
         depth_sc = as.numeric(scale(depth)),
         oxygen_sc = as.numeric(scale(oxygen)),
         density_saduria_sc = as.numeric(scale(density_saduria)),
         flounder_density_sc = as.numeric(scale(log(flounder_density))),
         large_cod_density_sc = as.numeric(scale(log(large_cod_density))),
         small_cod_density_sc = as.numeric(scale(log(small_cod_density)))) |> 
  # Scale length by group ..
  mutate(pred_length_cm_sc = as.numeric(scale(pred_length_cm)),
         .by = group)


d |> 
  rename("Relative benthic weight" = "benthic_rel_weight",
         "Relative Saduria weight" = "saduria_rel_weight") |> 
  pivot_longer(c("Relative benthic weight", "Relative Saduria weight")) |>
  mutate(group = str_to_sentence(group)) |> 
  ggplot(aes(value)) + 
  ggh4x::facet_grid2(factor(group, levels = c("Flounder", "Small cod", "Large cod"))~name,
                     scales = "free", independent = "all") + 
  geom_density(color = NA, alpha = 0.8, fill = "grey30") + 
  labs(y = "Density", x = "Value") + 
  NULL

ggsave(paste0(home, "/figures/supp/data_distribution.pdf"), width = 17, height = 17, units = "cm")

# Size distribution by year

t <- d |> filter(group == "large cod")

hist(t$pred_length_cm)

d |> 
  distinct(pred_id, .keep_all = TRUE) |> 
  mutate(group = str_to_sentence(group)) |> 
  ggplot(aes(pred_length_cm, color = as.factor(year))) + 
  facet_wrap(~group, scales = "free", ncol = 1) + 
  geom_density(alpha = 0.8, fill = NA) + 
  labs(y = "Density", x = "Predator length (cm)", color = "Year") + 
  scale_color_viridis(discrete = TRUE, option = "mako") +
  NULL

ggsave(paste0(home, "/figures/supp/predator_sizes.pdf"), width = 11, height = 17, units = "cm")

# Sample size per haul
d |> 
  summarise(n = length(unique(pred_id)), .by = haul_id) |> 
  ggplot(aes(n)) +
  geom_histogram()

# How many with 3 or fewer and what are their sample sizes?
d |> 
  summarise(n = length(unique(pred_id)), .by = haul_id) |> 
  mutate(s = ifelse(n < 3, "s", "l")) |> 
  summarise(n = n(), .by = s)

# A tibble: 2 × 2
  s         n
  <chr> <int>
1 l       316
2 s        29

# (29/316)*100

# What is the size range of those 10%?
d |> 
  mutate(n = length(unique(pred_id)), .by = haul_id) |> 
  filter(n > 1 & n < 11) |> 
  summarise(max = max(pred_length_cm),
            min = min(pred_length_cm),
            .by = c(group, haul_id)) |> 
  mutate(diff = max - min) |> 
  as.data.frame()

       group    haul_id max min diff
1  large cod   2015_4_6  49  46    3
2  large cod  2015_4_12  49  31   18
3  large cod  2015_4_32  38  27   11
4  small cod  2015_4_32  18  16    2
5  small cod  2015_4_51  21  21    0
6  large cod  2015_4_51  45  45    0
7  large cod  2015_4_53  45  41    4
8   flounder   2015_4_6  39  35    4
9   flounder  2015_4_12  26  22    4
10  flounder  2015_4_32  31   9   22
11 small cod  2016_1_16  18  16    2
12 large cod  2016_1_23  46  29   17
13 large cod  2016_1_37  49  30   19
14 small cod  2016_1_37  25  25    0
15 large cod  2016_1_53  34  31    3
16 large cod  2016_1_73  48  26   22
17 small cod  2016_1_73  22  18    4
18 large cod  2016_1_81  47  47    0
19 small cod  2016_1_81  14  14    0
20 small cod  2016_1_87  17  13    4
21 large cod  2016_1_89  50  48    2
22 small cod  2016_4_17  25  25    0
23 small cod  2016_4_62  14  12    2
24 large cod  2016_4_62  47  47    0
25  flounder  2016_4_17  29  29    0
26  flounder  2016_4_62  35  35    0
27 small cod  2017_1_17  12  11    1
28 large cod  2017_1_60  46  45    1
29 large cod  2017_1_95  45  45    0
30 small cod  2017_1_95  13  12    1
31  flounder  2017_1_95  36  15   21
32  flounder  2017_1_99  38  35    3
33  flounder  2017_1_40  34  20   14
34 small cod  2017_4_52  21  21    0
35 large cod  2017_4_52  28  28    0
36  flounder  2017_4_52  31  20   11
37 large cod  2018_1_52  44  26   18
38 small cod  2018_1_52  25  24    1
39 large cod  2019_4_91  36  29    7
40  flounder  2019_4_82  24  21    3
41  flounder  2019_4_91  26  21    5
42 large cod  2020_1_65  27  27    0
43  flounder  2020_1_57  23  15    8
44  flounder  2020_1_65  26  21    5
45 small cod 2020_4_213  24  24    0
46  flounder 2020_4_213  28  20    8
47  flounder 2020_4_214  27  22    5
48 large cod  2021_1_56  33  26    7
49 small cod  2021_1_56  23  20    3
50 large cod  2021_1_83  28  26    2
51 large cod  2021_1_84  31  31    0
52 small cod  2021_1_84  22  10   12
53 small cod  2021_1_85  12  10    2
54 large cod  2021_1_86  32  29    3
55 small cod  2021_1_86  10  10    0
56 large cod  2021_1_89  29  27    2
57 large cod  2021_1_91  37  37    0
58 small cod  2021_1_91  24  21    3
59 large cod  2021_1_98  40  27   13
60 small cod  2021_1_98  23  23    0
61 large cod  2021_1_99  40  27   13
62 small cod  2021_1_99  25  19    6
63 large cod 2021_1_100  38  38    0
64 small cod 2021_1_100  20  20    0
65 large cod 2021_1_103  44  38    6
66 small cod 2021_1_103  25  10   15
67 large cod 2021_1_107  36  33    3
68 large cod 2021_1_108  41  31   10
69 small cod 2021_1_108  10  10    0
70 large cod 2021_4_258  34  34    0
71 small cod 2021_4_258  23  23    0
72  flounder 2021_4_258  25  18    7
73  flounder  2022_1_88  26  22    4
74 large cod 2022_4_241  32  32    0
75 large cod 2022_4_266  40  29   11
76 small cod 2022_4_266  23  19    4
77 large cod 2022_4_275  34  29    5
78  flounder 2022_4_241  34  20   14
79  flounder 2022_4_266  24  23    1

  #summarise(mean = mean(diff))

Quick explore

Correlation between variables

# Plot correlation between variables
d_cor <- d |>
  dplyr::select("oxygen_sc", "density_saduria_sc", "flounder_density_sc",
                "large_cod_density_sc", "small_cod_density_sc", "depth_sc")

corr <- round(cor(d_cor), 1)

library(ggcorrplot)
ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.3))

Sample size

d |>
  group_by(species) |> 
  summarise(n = n())

group_by: one grouping variable (species)

summarise: now 2 rows and 2 columns, ungrouped

# A tibble: 2 × 2
  species      n
  <chr>    <int>
1 Cod       5484
2 Flounder  3851

d |>
  group_by(species, quarter) |> 
  summarise(n = n())

group_by: 2 grouping variables (species, quarter)

summarise: now 4 rows and 3 columns, one group variable remaining (species)

# A tibble: 4 × 3
# Groups:   species [2]
  species  quarter     n
  <chr>      <dbl> <int>
1 Cod            1  2906
2 Cod            4  2578
3 Flounder       1  2081
4 Flounder       4  1770

Fit models

Groups are: small cod, large cod and flounder. Response variables are: saduria_rel_weight, benthic_rel_weight or total weight. The latter is only for adult cod, because essentially all prey are benthic for small cod and flounder.

Model random effect structure is selected with AIC (see script 02-)

# This is the reason we don't do total weight for flounder and small cod
d |>
  filter(tot_rel_weight > 0) |> 
  group_by(group) |> 
  mutate(ben_prop = benthic_rel_weight / tot_rel_weight) |> 
  summarise(mean_ben_prop = mean(ben_prop))

# A tibble: 3 × 2
  group     mean_ben_prop
  <chr>             <dbl>
1 flounder          0.978
2 large cod         0.591
3 small cod         0.956

Covariates are: ~ 0 + fyear + fquarter + depth_sc + spatial + spatiotemporal random fields + density covariates. For saduria, we use saduria also in interaction with cod and flounder. For cod we use small cod because large and small cod are very correlated. For benthic and total prey, we instead use oxygen, more as a proxy, as the interaction variable

Main models

pred_flounder_sad <- list()
pred_flounder_ben <- list()
pred_cod_sad <- list()
pred_cod_ben <- list()
coef_sad <- list()
coef_ben <- list()
res_sad <- list()
res_ben <- list()
random_sad <- list()
random_ben <- list()
range_sad <- list()
range_ben <- list()

for(i in unique(d$group)) {
  
  dd <- filter(d, group == i)
  
  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 5)

  ggplot() +
    inlabru::gg(mesh$mesh) +
    coord_fixed() +
    geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
    annotate("text", -Inf, Inf, label = paste("n knots =", mesh$mesh$n), hjust = -0.3, vjust = 3) + 
    labs(x = "Easting (km)", y = "Northing (km)")
  
  ggsave(paste0(home, "/figures/supp/mesh_", i, ".pdf"), width = 17, height = 17, units = "cm")
  
  
  # Saduria model
  
  m_sad <- sdmTMB(saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + pred_length_cm_sc +
                    small_cod_density_sc*density_saduria_sc + 
                    flounder_density_sc*density_saduria_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(),
                  spatiotemporal = "iid", 
                  spatial = "on",
                  time = "year")
  print(i)
  sanity(m_sad)
  print(m_sad)
  
  
  # Benthic model
  
    if(unique(dd$group) %in% c("large cod", "small cod")) {
      
      
      m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc +
                        small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                      data = dd,
                      mesh = mesh,
                      family = tweedie(),
                      spatiotemporal = "iid", 
                      spatial = "off", 
                      time = "year")
      print(i)
      sanity(m_ben)
      print(m_ben)
      
    } else {
      
      m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc +
                        small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                      data = dd,
                      mesh = mesh,
                      family = tweedie(),
                      spatiotemporal = "iid", 
                      spatial = "on",
                      time = "year")
      print(i)
      sanity(m_ben)
      print(m_ben)
      
    }
       
   
  # Spatial and spatiotemporal random effects
  d_haul <- dd |>
    distinct(haul_id, .keep_all = TRUE)

  preds_sad <- predict(m_sad, newdata = d_haul)
  preds_ben <- predict(m_ben, newdata = d_haul)

  random_sad[[i]] <- preds_sad
  random_ben[[i]] <- preds_ben

  # Residuals
  samps <- sdmTMBextra::predict_mle_mcmc(m_sad, mcmc_iter = 401, mcmc_warmup = 400)
  mcmc_res <- residuals(m_sad, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)

  res_sad[[i]] <- dd

  samps <- sdmTMBextra::predict_mle_mcmc(m_ben, mcmc_iter = 401, mcmc_warmup = 400)
  mcmc_res <- residuals(m_ben, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)

  res_ben[[i]] <- dd


  # Ranges
  range_sad[[i]] <- tidy(m_sad, effects = "ran_pars") |> filter(term == "range") |> mutate(group = i, model = "saduria")
  range_ben[[i]] <- tidy(m_ben, effects = "ran_pars") |> filter(term == "range") |> mutate(group = i, model = "benthos")


  # Conditional effects: flounder
  nd_flounder <- data.frame(expand_grid(
    density_saduria_sc = c(quantile(d$density_saduria_sc, probs = 0.05),
                           mean(d$density_saduria_sc),
                           quantile(d$density_saduria_sc, probs = 0.95)),
    flounder_density_sc = seq(quantile(dd$flounder_density_sc, probs = 0.05),
                              quantile(dd$flounder_density_sc, probs = 0.95),
                              length.out = 50))) |>
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           pred_length_cm_sc = 0,
           oxygen_sc = 0,
           depth_sc = 0,
           small_cod_density_sc = 0)

  preds_flounder_sad <- predict(m_sad, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  preds_flounder_ben <- predict(m_ben, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  pred_flounder_sad[[i]] <- preds_flounder_sad |> mutate(group = i, xvar = "flounder")
  pred_flounder_ben[[i]] <- preds_flounder_ben |> mutate(group = i, xvar = "flounder")

  # Conditional effects: cod
  nd_cod <- data.frame(expand_grid(
    density_saduria_sc = c(quantile(d$density_saduria_sc, probs = 0.05),
                           quantile(d$density_saduria_sc, probs = 0.95)),
    small_cod_density_sc = seq(quantile(dd$small_cod_density_sc, probs = 0.05),
                               quantile(dd$small_cod_density_sc, probs = 0.95),
                               length.out = 50))) |>
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           pred_length_cm_sc = 0,
           oxygen_sc = 0,
           depth_sc = 0,
           flounder_density_sc = 0) #

  preds_cod_sad <- predict(m_sad, newdata = nd_cod, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  preds_cod_ben <- predict(m_ben, newdata = nd_cod, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  pred_cod_sad[[i]] <- preds_cod_sad |> mutate(group = i, xvar = "cod")
  pred_cod_ben[[i]] <- preds_cod_ben |> mutate(group = i, xvar = "cod")

  # Coefficients
  coefs_sad <- bind_rows(tidy(m_sad, effects = "fixed", conf.int = TRUE)) |>
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))

  coefs_ben <- bind_rows(tidy(m_ben, effects = "fixed", conf.int = TRUE)) |>
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))

  coef_sad[[i]] <- coefs_sad |> mutate(group = i)
  coef_ben[[i]] <- coefs_ben |> mutate(group = i)

}

filter: removed 5,801 rows (62%), 3,534 rows remaining

Loading required namespace: INLA

[1] "large cod"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     pred_length_cm_sc + small_cod_density_sc * density_saduria_sc + 
 Formula:     flounder_density_sc * density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                                        coef.est coef.se
fyear2015                                  -8.57    0.87
fyear2016                                 -10.06    0.71
fyear2017                                 -10.71    0.74
fyear2018                                 -10.08    0.81
fyear2019                                 -11.86    1.03
fyear2020                                 -10.48    0.70
fyear2021                                 -10.82    0.73
fyear2022                                 -11.37    0.81
fquarter4                                  -0.37    0.27
depth_sc                                   -0.83    0.27
oxygen_sc                                   0.19    0.22
pred_length_cm_sc                          -0.60    0.08
small_cod_density_sc                        0.28    0.20
density_saduria_sc                          0.35    0.24
flounder_density_sc                        -0.38    0.19
small_cod_density_sc:density_saduria_sc     0.18    0.20
density_saduria_sc:flounder_density_sc      0.25    0.17

Dispersion parameter: 0.11
Tweedie p: 1.45
Matérn range: 29.80
Spatial SD: 2.04
Spatiotemporal IID SD: 1.70
ML criterion at convergence: -807.559

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "large cod"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc + 
 Formula:     small_cod_density_sc * oxygen_sc + flounder_density_sc * 
 Formula:     oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                               coef.est coef.se
fyear2015                         -6.42    0.30
fyear2016                         -6.65    0.19
fyear2017                         -6.60    0.19
fyear2018                         -6.51    0.22
fyear2019                         -7.30    0.29
fyear2020                         -6.74    0.19
fyear2021                         -6.54    0.18
fyear2022                         -6.75    0.22
fquarter4                          0.61    0.11
depth_sc                          -0.25    0.07
pred_length_cm_sc                 -0.30    0.03
small_cod_density_sc              -0.06    0.06
oxygen_sc                          0.19    0.06
flounder_density_sc                0.04    0.06
small_cod_density_sc:oxygen_sc     0.00    0.03
oxygen_sc:flounder_density_sc      0.03    0.04

Dispersion parameter: 0.25
Tweedie p: 1.62
Matérn range: 9.59
Spatiotemporal IID SD: 1.37
ML criterion at convergence: -7817.046

See ?tidy.sdmTMB to extract these values as a data frame.

distinct: removed 3,230 rows (91%), 304 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.005223 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 52.23 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 57.706 seconds (Warm-up)
Chain 1:                0.097 seconds (Sampling)
Chain 1:                57.803 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00474 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 47.4 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 30.391 seconds (Warm-up)
Chain 1:                0.044 seconds (Sampling)
Chain 1:                30.435 seconds (Total)
Chain 1: 

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

filter: removed 3 rows (75%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'flounder_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA

filter: removed 7,385 rows (79%), 1,950 rows remaining

[1] "small cod"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     pred_length_cm_sc + small_cod_density_sc * density_saduria_sc + 
 Formula:     flounder_density_sc * density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                                        coef.est coef.se
fyear2015                                  -8.90    0.99
fyear2016                                  -9.43    0.86
fyear2017                                  -9.74    0.82
fyear2018                                  -9.67    1.05
fyear2019                                  -8.45    1.00
fyear2020                                  -9.34    0.80
fyear2021                                 -10.10    0.80
fyear2022                                  -9.65    0.84
fquarter4                                  -1.37    0.32
depth_sc                                   -1.05    0.34
oxygen_sc                                  -0.34    0.26
pred_length_cm_sc                           0.66    0.10
small_cod_density_sc                        0.46    0.37
density_saduria_sc                          0.56    0.31
flounder_density_sc                        -0.38    0.22
small_cod_density_sc:density_saduria_sc     0.72    0.36
density_saduria_sc:flounder_density_sc      0.23    0.19

Dispersion parameter: 0.17
Tweedie p: 1.50
Matérn range: 31.69
Spatial SD: 2.79
Spatiotemporal IID SD: 1.31
ML criterion at convergence: -640.781

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "small cod"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc + 
 Formula:     small_cod_density_sc * oxygen_sc + flounder_density_sc * 
 Formula:     oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                               coef.est coef.se
fyear2015                         -5.95    0.31
fyear2016                         -5.69    0.20
fyear2017                         -5.73    0.19
fyear2018                         -5.67    0.24
fyear2019                         -6.35    0.32
fyear2020                         -5.78    0.19
fyear2021                         -6.02    0.18
fyear2022                         -5.81    0.21
fquarter4                          0.48    0.10
depth_sc                          -0.35    0.07
pred_length_cm_sc                 -0.23    0.03
small_cod_density_sc               0.11    0.09
oxygen_sc                          0.16    0.06
flounder_density_sc                0.03    0.06
small_cod_density_sc:oxygen_sc    -0.05    0.07
oxygen_sc:flounder_density_sc      0.02    0.06

Dispersion parameter: 0.10
Tweedie p: 1.54
Matérn range: 19.53
Spatiotemporal IID SD: 0.84
ML criterion at convergence: -5599.649

See ?tidy.sdmTMB to extract these values as a data frame.

distinct: removed 1,669 rows (86%), 281 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.004663 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 46.63 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 61.381 seconds (Warm-up)
Chain 1:                0.084 seconds (Sampling)
Chain 1:                61.465 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.004289 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 42.89 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 29.327 seconds (Warm-up)
Chain 1:                0.077 seconds (Sampling)
Chain 1:                29.404 seconds (Total)
Chain 1: 

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

filter: removed 3 rows (75%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'flounder_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA

filter: removed 5,484 rows (59%), 3,851 rows remaining

[1] "flounder"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     pred_length_cm_sc + small_cod_density_sc * density_saduria_sc + 
 Formula:     flounder_density_sc * density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                                        coef.est coef.se
fyear2015                                  -8.67    1.12
fyear2016                                  -6.23    0.94
fyear2017                                  -8.47    0.94
fyear2018                                  -8.75    1.00
fyear2019                                  -8.31    1.05
fyear2020                                  -8.42    0.95
fyear2021                                 -11.31    1.25
fyear2022                                  -9.10    0.99
fquarter4                                  -0.64    0.21
depth_sc                                   -0.37    0.22
oxygen_sc                                  -0.03    0.21
pred_length_cm_sc                           0.26    0.05
small_cod_density_sc                        0.45    0.13
density_saduria_sc                          0.26    0.16
flounder_density_sc                        -0.28    0.15
small_cod_density_sc:density_saduria_sc     0.09    0.10
density_saduria_sc:flounder_density_sc     -0.25    0.14

Dispersion parameter: 0.17
Tweedie p: 1.49
Matérn range: 72.51
Spatial SD: 2.17
Spatiotemporal IID SD: 1.59
ML criterion at convergence: -1808.407

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "flounder"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc + 
 Formula:     small_cod_density_sc * oxygen_sc + flounder_density_sc * 
 Formula:     oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                               coef.est coef.se
fyear2015                         -5.12    0.31
fyear2016                         -4.94    0.24
fyear2017                         -4.98    0.20
fyear2018                         -5.51    0.24
fyear2019                         -5.17    0.26
fyear2020                         -4.70    0.19
fyear2021                         -5.34    0.25
fyear2022                         -5.03    0.21
fquarter4                          0.04    0.10
depth_sc                          -0.37    0.10
pred_length_cm_sc                  0.01    0.02
small_cod_density_sc               0.22    0.06
oxygen_sc                         -0.05    0.07
flounder_density_sc               -0.03    0.08
small_cod_density_sc:oxygen_sc     0.01    0.03
oxygen_sc:flounder_density_sc     -0.03    0.05

Dispersion parameter: 0.17
Tweedie p: 1.50
Matérn range: 18.51
Spatial SD: 0.95
Spatiotemporal IID SD: 0.67
ML criterion at convergence: -5534.967

See ?tidy.sdmTMB to extract these values as a data frame.

distinct: removed 3,590 rows (93%), 261 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.005991 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 59.91 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 314.008 seconds (Warm-up)
Chain 1:                0.856 seconds (Sampling)
Chain 1:                314.864 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.005855 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 58.55 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 50.657 seconds (Warm-up)
Chain 1:                0.113 seconds (Sampling)
Chain 1:                50.77 seconds (Total)
Chain 1: 

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'flounder_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA

Now do a separate model for adult cod looking at total prey

  dd <- filter(d, group == "large cod")

filter: removed 5,801 rows (62%), 3,534 rows remaining

  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 5)

  # Total model
  # NOTE: turning off spatial here due to convergence and AIC
  m_tot <- sdmTMB(tot_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc +
                    large_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(),
                  spatiotemporal = "iid", 
                  spatial = "on",
                  time = "year")
  
  sanity(m_tot)

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

  print(m_tot)

Spatiotemporal model fit by ML ['sdmTMB']
Formula: tot_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc + 
 Formula:     large_cod_density_sc * oxygen_sc + flounder_density_sc * 
 Formula:     oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                               coef.est coef.se
fyear2015                         -4.30    0.23
fyear2016                         -4.58    0.16
fyear2017                         -4.66    0.16
fyear2018                         -4.41    0.18
fyear2019                         -4.85    0.23
fyear2020                         -4.58    0.15
fyear2021                         -4.33    0.15
fyear2022                         -4.62    0.17
fquarter4                         -0.09    0.09
depth_sc                          -0.04    0.06
pred_length_cm_sc                  0.27    0.03
large_cod_density_sc               0.11    0.07
oxygen_sc                         -0.06    0.05
flounder_density_sc                0.03    0.06
large_cod_density_sc:oxygen_sc     0.17    0.05
oxygen_sc:flounder_density_sc     -0.11    0.04

Dispersion parameter: 0.57
Tweedie p: 1.71
Matérn range: 25.60
Spatial SD: 0.29
Spatiotemporal IID SD: 0.51
ML criterion at convergence: -7957.724

See ?tidy.sdmTMB to extract these values as a data frame.

  # Residuals
  samps <- sdmTMBextra::predict_mle_mcmc(m_tot, mcmc_iter = 401, mcmc_warmup = 400)

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00597 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 59.7 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 49.144 seconds (Warm-up)
Chain 1:                0.103 seconds (Sampling)
Chain 1:                49.247 seconds (Total)
Chain 1: 

  mcmc_res <- residuals(m_tot, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)
    
  res_tot <- dd
  
  # Range
  range_tot <- tidy(m_tot, effects = "ran_pars") |> filter(term == "range") |> mutate(group = "large cod", model = "total")

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

  # Spatial and spatiotemporal random effects
  d_haul <- dd |> 
    distinct(haul_id, .keep_all = TRUE)

distinct: removed 3,230 rows (91%), 304 rows remaining

  preds_tot <- predict(m_tot, newdata = d_haul)
  
  # Coefficients
  coefs_tot <- bind_rows(tidy(m_tot, effects = "fixed", conf.int = TRUE)) |> 
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

  coefs_tot <- coefs_tot |> mutate(group = "large cod")

mutate: new variable 'group' (character) with one unique value and 0% NA

  # Conditional effects: flounder
  nd_flounder <- data.frame(expand_grid(
    flounder_density_sc = seq(quantile(dd$flounder_density_sc, probs = 0.05),
                              quantile(dd$flounder_density_sc, probs = 0.95),
                              length.out = 50))) |>
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           pred_length_cm_sc = 0,
           oxygen_sc = 0,
           depth_sc = 0,
           large_cod_density_sc = 0)

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'large_cod_density_sc' (double) with one unique value and 0% NA

  preds_flounder_tot <- predict(m_tot, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  pred_flounder_tot <- preds_flounder_tot |> mutate(group = "large cod", xvar = "flounder")

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

Make dataframes

coef_df <- bind_rows(bind_rows(coef_sad) |> mutate(model = "Saduria"),
                     bind_rows(coef_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

coef_df <- coef_df |> bind_rows(coefs_tot |> mutate(model = "Total"))

mutate: new variable 'model' (character) with one unique value and 0% NA

pred_cod_df <- bind_rows(bind_rows(pred_cod_sad) |> mutate(model = "Saduria"),
                         bind_rows(pred_cod_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

pred_flounder_df <- bind_rows(bind_rows(pred_flounder_sad) |> mutate(model = "Saduria"),
                              bind_rows(pred_flounder_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

res_df <- bind_rows(bind_rows(res_sad) |> mutate(model = "Saduria"),
                    bind_rows(res_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

res_df <- res_df |> bind_rows(res_tot |> mutate(model = "Total"))

mutate: new variable 'model' (character) with one unique value and 0% NA

random_df <- bind_rows(bind_rows(random_sad) |> mutate(model = "Saduria"),
                       bind_rows(random_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

random_df <- random_df |> bind_rows(preds_tot |> mutate(model = "Total"))

mutate: new variable 'model' (character) with one unique value and 0% NA

range_df <- bind_rows(range_tot, bind_rows(range_ben), bind_rows(range_sad))

Plot spatial random effects

random_df <- random_df |>
  mutate(group = str_to_sentence(group))

mutate: changed 1,996 values (100%) of 'group' (0 new NA)

# Saduria
plot_map_fc +
  geom_point(data = random_df |> filter(model == "Saduria"),
             aes(X*1000, Y*1000, color = omega_s), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~factor(group, levels = c("Flounder", "Small cod", "Large cod"))) +
  labs(color = "Spatial\nrandom effect") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        #legend.direction = "vertical",
        legend.key.width = unit(0.6, "cm"),
        legend.key.height = unit(0.2, "cm"))

filter: removed 1,150 rows (58%), 846 rows remaining

ggsave(paste0(home, "/figures/supp/omega_sad.pdf"), width = 17, height = 9, units = "cm")


# Now do benthos (only for flounder)
plot_map_fc +
  geom_point(data = random_df |> filter(model == "Benthos" & group == "Flounder"),
             aes(X*1000, Y*1000, color = omega_s), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ group) +
  labs(color = "Spatial\nrandom effect") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "right",
        legend.direction = "vertical",
        legend.key.width = unit(0.4, "cm"),
        legend.key.height = unit(0.4, "cm"))

filter: removed 1,735 rows (87%), 261 rows remaining

ggsave(paste0(home, "/figures/supp/omega_ben.pdf"), width = 11, height = 11, units = "cm")

Plot spatiotemporal random effects

# Saduria
sad_eps_sc <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Saduria" & group == "Small cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,715 rows (86%), 281 rows remaining

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

sad_eps_sc

ggsave(paste0(home, "/figures/supp/epsilon_sad_small_cod.pdf"), width = 17, height = 17, units = "cm")

sad_eps_lc <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Saduria" & group == "Large cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,692 rows (85%), 304 rows remaining

sad_eps_lc

ggsave(paste0(home, "/figures/supp/epsilon_sad_large_cod.pdf"), width = 17, height = 17, units = "cm")

sad_eps_f <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Saduria" & group == "Flounder"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,735 rows (87%), 261 rows remaining

sad_eps_f

ggsave(paste0(home, "/figures/supp/epsilon_sad_flounder.pdf"), width = 17, height = 17, units = "cm")


# Benthos
ben_eps_sc <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Benthos" & group == "Small cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,715 rows (86%), 281 rows remaining

ben_eps_sc

ggsave(paste0(home, "/figures/supp/epsilon_ben_small_cod.pdf"), width = 17, height = 17, units = "cm")

ben_eps_lc <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Benthos" & group == "Large cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,692 rows (85%), 304 rows remaining

ben_eps_lc

ggsave(paste0(home, "/figures/supp/epsilon_ben_large_cod.pdf"), width = 17, height = 17, units = "cm")

ben_eps_f <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Benthos" & group == "Flounder"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,735 rows (87%), 261 rows remaining

ben_eps_f

ggsave(paste0(home, "/figures/supp/epsilon_ben_flounder.pdf"), width = 17, height = 17, units = "cm")


# Total
tot_eps <- plot_map_fc +
    geom_point(data = preds_tot, aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
    scale_color_gradient2() +
    facet_wrap(~ year) +
    labs(color = "Spatiotemporal\nrandom effect") +
    theme(legend.position = c(0.84, 0.16),
          axis.text.x = element_text(angle = 90))

tot_eps

Plot range

#pal <- brewer.pal(n = 8, name = "Dark2")[c(2, 7, 6)]
pal <- (brewer.pal(n = 11, name = "RdYlBu")[c(11, 4, 1)])

range_df |> arrange(estimate)

# A tibble: 7 × 7
  term  estimate std.error conf.low conf.high group     model  
  <chr>    <dbl>     <dbl>    <dbl>     <dbl> <chr>     <chr>  
1 range     9.59      2.52     5.74      16.1 large cod benthos
2 range    18.5       4.39    11.6       29.5 flounder  benthos
3 range    19.5       5.59    11.2       34.2 small cod benthos
4 range    25.6       6.55    15.5       42.3 large cod total  
5 range    29.8      10.9     14.5       61.2 large cod saduria
6 range    31.7      10.3     16.8       59.8 small cod saduria
7 range    72.5      16.4     46.5      113.  flounder  saduria

range_df |> 
  mutate(group = str_to_sentence(group),
         model2 = ifelse(model == "benthos", "Benthic prey", model),
         model2 = ifelse(model == "saduria", "Saduria", model2),
         model2 = ifelse(model == "total", "Total prey", model2)) |> 
  ggplot(aes(model2, estimate, color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) + 
  geom_point(size = 2) +
  geom_hline(yintercept = 5, linetype = 2, alpha = 0.5) +
  scale_color_manual(values = pal) + 
  labs(x = "", y = "Range (km)", color = "Predator") + 
  theme(aspect.ratio = 1,
        legend.position = c(0.86, 0.86)) 

mutate: changed 7 values (100%) of 'group' (0 new NA)

        new variable 'model2' (character) with 3 unique values and 0% NA

ggsave(paste0(home, "/figures/supp/ranges.pdf"), width = 11, height = 11, units = "cm")

Plot residuals

# Plot residuals
res_df |> 
  mutate(group = str_to_title(group)) |> 
  ggplot(aes(sample = res)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_grid(model ~ group) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  theme(aspect.ratio = 1)

mutate: changed 22,204 values (100%) of 'group' (0 new NA)

ggsave(paste0(home, "/figures/supp/qq_relative_prey_weight.pdf"), width = 17, height = 17, units = "cm")

Plot coefficients

coef_df$term |> unique()

 [1] "fyear2015"                              
 [2] "fyear2016"                              
 [3] "fyear2017"                              
 [4] "fyear2018"                              
 [5] "fyear2019"                              
 [6] "fyear2020"                              
 [7] "fyear2021"                              
 [8] "fyear2022"                              
 [9] "fquarter4"                              
[10] "depth_sc"                               
[11] "oxygen_sc"                              
[12] "pred_length_cm_sc"                      
[13] "small_cod_density_sc"                   
[14] "density_saduria_sc"                     
[15] "flounder_density_sc"                    
[16] "small_cod_density_sc:density_saduria_sc"
[17] "density_saduria_sc:flounder_density_sc" 
[18] "small_cod_density_sc:oxygen_sc"         
[19] "oxygen_sc:flounder_density_sc"          
[20] "large_cod_density_sc"                   
[21] "large_cod_density_sc:oxygen_sc"         

# Fix some names
coef_df2 <- coef_df |>
  filter(!grepl('year', term)) |>
  filter(!grepl('quarter', term)) |>
  mutate(term = str_remove_all(term, "_sc"),
         term = str_remove_all(term, "density"),
         term = str_replace_all(term, "_", ""),
         term = str_replace_all(term, "geco", "ge co"),
         term = str_replace_all(term, "llco", "ll co"),
         term = str_replace(term, ":", " × "),
         term = str_to_sentence(term),
         group = str_to_sentence(group)) |> 
  mutate(model2 = ifelse(model == "Saduria", "Prey=Saduria", NA),
         model2 = ifelse(model == "Benthos", "Prey=Benthic prey", model2),
         model2 = ifelse(model == "Total", "Prey=Total prey", model2)) |> 
  mutate(sig2 = ifelse(sig == "Y", "N", "Y")) |> 
  mutate(term = ifelse(term == "Predlengthcm", "Predator length", term))

filter: removed 56 rows (49%), 59 rows remaining

filter: removed 7 rows (12%), 52 rows remaining

mutate: changed 52 values (100%) of 'term' (0 new NA)

        changed 52 values (100%) of 'group' (0 new NA)

mutate: new variable 'model2' (character) with 3 unique values and 0% NA

mutate: new variable 'sig2' (character) with 2 unique values and 0% NA

mutate: changed 7 values (13%) of 'term' (0 new NA)

p1 <- 
  coef_df2 |> 
  filter(model == "Benthos") |> 
  ggplot(aes(estimate, term, alpha = sig2,
                 color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
  geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
  facet_wrap2(~model2, ncol = 2, scales = "free") +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(1, 0.4)) +
  scale_color_manual(values = pal) +
  labs(x = "", y = "", alpha = "95% CI crossing 0", color = "Predator") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
         alpha = guide_legend(title.position = "top", title.hjust = 0.5)) +
  theme(legend.position = c(0.75, 0.2),
        legend.direction = "vertical",
        legend.box.spacing = unit(-3, "pt"),
        legend.margin = margin(0, 0, 0, 0)) +
  NULL

filter: removed 31 rows (60%), 21 rows remaining

p2 <- 
  coef_df2 |> 
  filter(model == "Saduria") |> 
  ggplot(aes(estimate, term, alpha = sig2,
                 color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
  geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
  facet_wrap2(~model2, ncol = 2, scales = "free") +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(1, 0.4)) +
  scale_color_manual(values = pal) +
  labs(x = "", y = "") +
  guides(color = "none",
         alpha = "none")

filter: removed 28 rows (54%), 24 rows remaining

p3 <- 
  coef_df2 |> 
  filter(model == "Total") |> 
  ggplot(aes(estimate, term, alpha = sig2,
                 color = "Large cod")) +
  geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
  facet_wrap2(~model2, ncol = 2, scales = "free") +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(1, 0.4)) +
  scale_color_manual(values = pal[3]) +
  labs(x = "", y = "") +
  guides(color = "none",
         alpha = "none")

filter: removed 45 rows (87%), 7 rows remaining

p2 / p1 / p3 +
  plot_annotation(tag_levels = "a") +
  plot_layout(axes = "collect", guides = "collect") & 
  theme(legend.position = "bottom") &
  coord_cartesian(xlim = c(-1.55, 1.15))

ggsave(paste0(home, "/figures/coefs.pdf"), width = 11, height = 22, units = "cm")

Plot year and quarter coefficients

# Fix some names
coef_df3 <- coef_df |>
  filter(grepl('year', term)) |>
  mutate(term = str_remove_all(term, "fyear"),
         group = str_to_sentence(group),
         term = as.numeric(term),
         model2 = ifelse(model == "Benthos", "Benthic prey", model),
         model2 = ifelse(model == "Saduria", "Saduria", model2),
         model2 = ifelse(model == "Total", "Total prey", model2))

filter: removed 59 rows (51%), 56 rows remaining

mutate: converted 'term' from character to double (0 new NA)

        changed 56 values (100%) of 'group' (0 new NA)

        new variable 'model2' (character) with 3 unique values and 0% NA

ggplot(coef_df3, aes(term, exp(estimate), color = factor(group, levels = c("Flounder", "Small cod", "Large cod")), 
                     fill = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
  facet_wrap(~model2, scales = "free", ncol = 1) +
  #geom_line(position = position_dodge(width = 0.5)) +
  geom_point(position = position_dodge(width = 0.4)) +
  #geom_ribbon(aes(ymin = exp(conf.low), ymax = exp(conf.high)), alpha = 0.2, color = NA) +
  geom_errorbar(aes(ymin = exp(conf.low), ymax = exp(conf.high)), alpha = 0.4, width = 0,
                position = position_dodge(width = 0.4)) +
  scale_color_manual(values = pal) +
  scale_fill_manual(values = pal) +
  labs(x = "Year", y = "Standardized coefficient", color = "") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5, ncol = 3),
         fill = "none") +
  theme(legend.position = c(0.5, 0.99),
        legend.direction = "vertical",
        legend.box.spacing = unit(-3, "pt"),
        legend.margin = margin(0, 0, 0, 0),
        strip.text.x.top = element_text(angle = 0, hjust = 0)) +
  NULL

ggsave(paste0(home, "/figures/supp/coefs_year.pdf"), width = 11, height = 21, units = "cm")

# Now do quarter
coef_df5 <- coef_df |>
  filter(term %in% c("fquarter4")) |> 
  mutate(group = str_to_sentence(group),
         model2 = ifelse(model == "Benthos", "Benthic prey", model),
         model2 = ifelse(model == "Saduria", "Saduria", model2),
         model2 = ifelse(model == "Total", "Total prey", model2)) |> 
  mutate(sig2 = ifelse(sig == "Y", "N", "Y"))

filter: removed 108 rows (94%), 7 rows remaining

mutate: changed 7 values (100%) of 'group' (0 new NA)

        new variable 'model2' (character) with 3 unique values and 0% NA

mutate: new variable 'sig2' (character) with 2 unique values and 0% NA

ggplot(coef_df5, aes(estimate, model2,
                     alpha = sig2,
                     color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(1, 0.4)) +
  scale_color_manual(values = pal) +
  labs(x = "", y = "Quarter 4 effect", alpha = "95% CI crossing 0", color = "Predator") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
         alpha = guide_legend(title.position = "top", title.hjust = 0.5)) +
  geom_stripped_rows(aes(y = model2), inherit.aes = FALSE) +
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box = "horizontal",
        legend.box.spacing = unit(-3, "pt"),
        legend.margin = margin(0, 0, 0, 0)) + 
  coord_cartesian(expand = 0)

ggsave(paste0(home, "/figures/supp/coefs_quarter.pdf"), width = 17, height = 11, units = "cm")

Conditional effects

# Which CI?
# https://www.calculator.net/confidence-interval-calculator.html
pred_df <- bind_rows(pred_cod_df, pred_flounder_df) |>
  mutate(group = str_to_sentence(group),
         sad = ifelse(density_saduria_sc == min(density_saduria_sc), "Low", NA),
         sad = ifelse(density_saduria_sc == max(density_saduria_sc), "High", sad)) |> 
  drop_na(sad)

mutate: changed 1,500 values (100%) of 'group' (0 new NA)

        new variable 'sad' (character) with 3 unique values and 20% NA

drop_na: removed 300 rows (20%), 1,200 rows remaining

pred_df2 <- bind_rows(pred_flounder_df,
                      pred_flounder_tot |> mutate(model = "Total")) |> 
  mutate(group = str_to_sentence(group),
         density_saduria_sc = replace_na(density_saduria_sc,
                                         median(density_saduria_sc, na.rm = TRUE)))

mutate: new variable 'model' (character) with one unique value and 0% NA

mutate: changed 50 values (5%) of 'density_saduria_sc' (50 fewer NA)

        changed 950 values (100%) of 'group' (0 new NA)

# 75% CI!!
ggplot(pred_df |> filter(model == "Saduria" & xvar == "flounder"),# & !group == "Large cod"),
       aes(flounder_density_sc, exp(est), color = sad, fill = sad)) +
  geom_ribbon(aes(ymin = exp(est - 1.150*est_se), ymax = exp(est + 1.150*est_se)),
              alpha = 0.3, color = NA) +
  geom_line() +
  facet_wrap(~factor(group, levels = c("Flounder", "Small cod", "Large cod")),
             scales = "free", 
             ncol = 3) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  coord_cartesian(#xlim = c(-0.5, 0.5), 
                  expand = 0, clip = "off") +
  labs(x = "Flounder density", y = "Relative saduria weight",
       color = "Saduria", fill = "Saduria") + 
  theme(legend.position = c(0.95, 0.84),
        strip.text.x.top = element_text(angle = 0, hjust = 0)) + 
  scale_x_continuous(breaks = c(-1, 0, 1)) +
  NULL

filter: removed 900 rows (75%), 300 rows remaining

ggsave(paste0(home, "/figures/conditional_saduria_flounder.pdf"), width = 19, height = 7, units = "cm")

Showing conditional effects of oxygen on small cod feeding on benthos

  dd <- filter(d, group == "small cod")

filter: removed 7,385 rows (79%), 1,950 rows remaining

  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 5)

  # Benthic model
  m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + 
                    small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(),
                  spatiotemporal = "iid",
                  spatial = "off",
                  time = "year")

  sanity(m_ben)

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

  nd <- data.frame(oxygen = seq(quantile(d$oxygen, probs = 0.05), quantile(d$oxygen, probs = 0.95),
                                length.out = 50)) |>
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           density_saduria_sc = 0,
           flounder_density_sc = 0,
           depth_sc = 0,
           small_cod_density_sc = 0) |>
    mutate(oxygen_sc = (oxygen - mean(d$oxygen)) / sd(d$oxygen))

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'density_saduria_sc' (double) with one unique value and 0% NA

        new variable 'flounder_density_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'oxygen_sc' (double) with 50 unique values and 0% NA

  p <- predict(m_ben, newdata = nd, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  ggplot(p, aes(oxygen, exp(est))) +
    geom_line() +
    theme_sleek(base_size = 14) +
    geom_hline(yintercept = 0.0025, col = "red") +
    geom_hline(yintercept = 0.0040, col = "red") +
    geom_vline(xintercept = 4.8, col = "red") +
    geom_vline(xintercept = 7.6, col = "red") +
    NULL